Read and preprocess the source data

Read raw discharge data

The starting point for hydrograph analysis is to obtain the source data. Let’s see how it goes with sample spas-zagorye.txt discharge data for \(1956-2020\) year range provided with grwat. This data is for Spas-Zagorye gauge on Protva river in Central European plane:

library(sf) # reading and manipulating spatial data
library(tidyverse) # general data wrangling
library(mapview) # interactive mapping of spatial data
library(ecmwfr) # this is to access ERA5 reanalysis data
library(grwat)

mapviewOptions(fgb = FALSE)

# this is path to sample data installed with grwat
path = system.file("extdata", "spas-zagorye.txt", package = "grwat")

# for your own data just provide the full path:
# path = /this/is/the/path/to/discharge/discharge.csv

hdata = read_delim(path, col_names = c('d', 'm', 'y', 'q'), col_types = 'iiid', delim = ' ') # read gauge data
head(hdata) # see the data
#> # A tibble: 6 x 4
#>       d     m     y     q
#>   <int> <int> <int> <dbl>
#> 1     1     1  1956  5.18
#> 2     2     1  1956  5.18
#> 3     3     1  1956  5.44
#> 4     4     1  1956  5.44
#> 5     5     1  1956  5.44
#> 6     6     1  1956  5.58

This example contains the minimum amount of data needed to use the grwat. The first three columns encode the date, and the last column is the main variable — discharge. If the date is encoded as a single column instead of three separate columns, this is also perfect for grwat as input:

hdata_alt = hdata %>% 
  transmute(D = lubridate::make_date(y, m, d), 
            Q = q)
head(hdata_alt)
#> # A tibble: 6 x 2
#>   D              Q
#>   <date>     <dbl>
#> 1 1956-01-01  5.18
#> 2 1956-01-02  5.18
#> 3 1956-01-03  5.44
#> 4 1956-01-04  5.44
#> 5 1956-01-05  5.44
#> 6 1956-01-06  5.58

Join meteorological variables

A sole discharge data is enough to separate the hydrograph into quickflow and baseflow, but is not sufficient to predict the genesis of quickflow cases. Was it due to rain or thaw? To answer such questions you also need precipitation and temperature data. Ideally, these must be measured at the gauge. But often such data is not available. In this case you need to mine this data from external sources.

One of the ways to obtain the temperature and precipitation data is to use reanalyses such as ERA5. Reanalysis data is arranged as regular grids with specific resolution. In particular, the ERA5 data has \(31\) km or \(0.28125\) degrees resolution. To use such data you must tolerate the fact that none of the reanalysis grid nodes will coincide with your gauge. Instead, you have to use the data, which is either

the closest reanalysis point
all points falling within the buffer zone of the gauge
all points falling within the basin of the the gauge.
all points falling within the buffer zone of the basin of the the gauge.

The last two options are the most adequate. Let’s see how it can be done.

Basin, buffered basin or buffered gauge

First, we need to read the basin spatial data:

# this is path to sample basin geopackage installed with grwat
path = system.file("extdata", "spas-zagorye.gpkg", package = "grwat")

# for your own data just provide the full path:
# path = /this/is/the/path/to/discharge/basin.shp

basin = st_read(path, layer = 'basin') # read basin region
#> Reading layer `basin' from data source 
#>   `/private/var/folders/5s/rkxr4m8j24569d_p6nj9ld200000gn/T/RtmpFIW0dE/temp_libpathe4ef26265f1d/grwat/extdata/spas-zagorye.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 1 feature and 7 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 35.41204 ymin: 54.88195 xmax: 36.84138 ymax: 55.57005
#> Geodetic CRS:  WGS 84
gauge = st_read(path, layer = 'gauge') # read gauge point
#> Reading layer `gauge' from data source 
#>   `/private/var/folders/5s/rkxr4m8j24569d_p6nj9ld200000gn/T/RtmpFIW0dE/temp_libpathe4ef26265f1d/grwat/extdata/spas-zagorye.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 36.62272 ymin: 55.03669 xmax: 36.62272 ymax: 55.03669
#> Geodetic CRS:  WGS 84
mapview(basin) + mapview(gauge)

Next, we can buffer the data on the specified distance to catch more reanalysis data. For this we use gr_buffer_geo function, which approximates geographic buffer of the specified radius:

basin_buffer = gr_buffer_geo(basin, 25000) 
mapview(basin_buffer, col.regions = 'red') +
  mapview(basin)

Alternatively you can just buffer the gauge point, though it is less meaningful since you will grab reanalysis data that falls out of gauge’s basin:

gauge_buffer = gr_buffer_geo(gauge, 50000) 
mapview(gauge_buffer, col.regions = 'red') +
  mapview(gauge)

Joining the East-European Plane reanalysis

grwat is packaged with daily reanalysis which covers the East European territory of Russia. It has a spatial resolution of \(0.75^{\circ} \times 0.75^{\circ}\), and temporal resolution of \(1\) (one) day. Data sources include:

CIRES-DOE from 01-01-1880 till 31-12-1978
ERA5 from 01-01-1979 till now.

The coverage of the reanalysis is shown below:

The reanalysis consists of two data files, each file is about 850 Mb in size:

Download this data using the FTP link for using with grwat.

rean = gr_read_rean('/Volumes/Data/Spatial/Reanalysis/grwat/pre_1880-2021.nc',
                    '/Volumes/Data/Spatial/Reanalysis/grwat/temp_1880-2021.nc') # read reanalysis data
hdata_rean = gr_join_rean(hdata, rean, basin_buffer) # join reanalysis data to hydrological series
#> Joining, by = c("Year", "Month", "Day")

head(hdata_rean$df)
#>   Day Month Year    Q   Temp  Prec
#> 1   1     1 1956 5.18  -6.46 0.453
#> 2   2     1 1956 5.18 -11.41 0.825
#> 3   3     1 1956 5.44 -10.74 0.260
#> 4   4     1 1956 5.44  -8.05 0.397
#> 5   5     1 1956 5.44 -11.73 0.102
#> 6   6     1 1956 5.58 -20.13 0.032

After reanalysis data are joined you can easily plot a map of the derived spatial configuration:

 # plot spatial configuration
m = mapview(basin_buffer, col.regions = 'red') +
  mapview(basin) +
  mapview(hdata_rean$pts, col.regions = 'black') +
  mapview(rean$pts, cex = 1)

box = st_bbox(basin_buffer)
center = st_coordinates(st_centroid(basin_buffer))
#> Warning in st_centroid.sf(basin_buffer): st_centroid assumes attributes are
#> constant over geometries of x

m@map %>% leaflet::setView(center[1], center[2], zoom = 7)

Joining the ERA5 reanalysis

To be described

Fill gaps

gr_fill_gaps() function allows to fill the gaps in the data using simple linear interpolation between nearest observations. You can limit the maximum gap extent by threshold autocorrelation value (autocorr parameter) or expliсit number of observations (nobserv parameter):

tab = gr_fill_gaps(hdata_rean$df,
                   nobserv = 10)
#> grwat: filled 0 observations using 10 days window

tab = gr_fill_gaps(hdata_rean$df,
                   autocorr = 0.7)
#> grwat: filled 0 observations using 5 days window

How the time window is computed in the second case? When the autocorr parameter is used, the ACF (autocorrelation function) is computed first, and then its values are used to obtain the time shift (in days) during which the autocorrelation is higher or equal to the specified value. For Spas-Zagorye data the time lag for \(ACF \geq 0.7\) is 5 days, whoch can be extracted using the purrr::detect_index() function:

afun = acf(hdata_rean$df$Q)
abline(h = 0.7, col = 'red')
abline(v = purrr::detect_index(afun$acf, ~ .x < 0.7), col = 'blue', lwd = 2)

Separate and summarize

Basic separation: quickflow and baseflow

The basic separation procedure is provided by get_baseflow() function. One of the most commonly used approaches is the method by Lyne-Hollick (1979):

resbase = tab %>% 
  mutate(Date = lubridate::dmy(paste(Day, Month, Year))) %>% 
  mutate(Qbase = gr_baseflow(Q, method = 'lynehollick'))

# quick look at the table
head(resbase, 10)
#> # A tibble: 10 x 8
#>      Day Month  Year     Q  Temp  Prec Date       Qbase
#>    <int> <dbl> <dbl> <dbl> <dbl> <dbl> <date>     <dbl>
#>  1     1     1  1956   5.2  -6.5   0.5 1956-01-01  3.70
#>  2     2     1  1956   5.2 -11.4   0.8 1956-01-02  3.79
#>  3     3     1  1956   5.4 -10.7   0.3 1956-01-03  3.88
#>  4     4     1  1956   5.4  -8.1   0.4 1956-01-04  3.96
#>  5     5     1  1956   5.4 -11.7   0.1 1956-01-05  4.04
#>  6     6     1  1956   5.6 -20.1   0   1956-01-06  4.12
#>  7     7     1  1956   5.6 -16.6   0   1956-01-07  4.19
#>  8     8     1  1956   5.4 -15.2   0   1956-01-08  4.26
#>  9     9     1  1956   5.4 -15.9   0.1 1956-01-09  4.32
#> 10    10     1  1956   5.3 -13.3   0.1 1956-01-10  4.39
  
resbase %>% 
  filter(Year == 2020) %>% 
  ggplot() +
  geom_area(aes(Date, Q), fill = 'steelblue') +
  geom_area(aes(Date, Qbase), fill = 'orangered')

The separation is mainly parameterized by smoothing parameter which is a = 0.925 by default, and number of passes, which are passes = 3 by default. Changing them affects the shape of a baseflow component:

resbase = tab %>% 
  mutate(Date = lubridate::dmy(paste(Day, Month, Year))) %>% 
  mutate(Qbase = gr_baseflow(Q, method = 'lynehollick', a = 0.95, passes = 5))
  
resbase %>% 
  filter(Year == 2020) %>% 
  ggplot() +
  geom_area(aes(Date, Q), fill = 'steelblue') +
  geom_area(aes(Date, Qbase), fill = 'orangered')

Let’s see how different separation methods act in comparison to each other:

methods = c("maxwell",
            "boughton",
            "jakeman",
            "lynehollick",
            "chapman")

plots = lapply(methods, function(m) {
  resbase = tab %>% 
    mutate(Date = lubridate::dmy(paste(Day, Month, Year))) %>% 
    mutate(Qbase = gr_baseflow(Q, method = m))

  resbase %>%
    filter(Year == 2020) %>%
    ggplot() +
    geom_area(aes(Date, Q), fill = 'steelblue') +
    geom_area(aes(Date, Qbase), fill = 'orangered') +
    labs(title = m)
})

patchwork::wrap_plots(plots, ncol = 2)

Advanced separation: genetic components

Advanced separation aims at revealing the genesis of the quickflow. Was it due to a rain or snowmelting? For this, a joint analysis of discharge, temperature and precipitation time series is performed by specialized algorithm, available from gr_separate() function in grwat package.

The genetic separation of hydrograph is controlled by more than 30 parameters. The names and meaning of these paramters can be learned thanks to gr_help_params() function:

gr_help_params()
#> # A tibble: 31 x 11
#>    number name   example desc    units formula comments   pics  desc_rus  role_1
#>     <dbl> <chr>  <chr>   <chr>   <chr> <chr>   <chr>      <lgl> <chr>     <chr> 
#>  1      1 mome   11      the fi… <NA>   <NA>   "not used… NA    месяц, с… <NA>  
#>  2      2 grad   1.5     rate o… %/day  <NA>   "if chang… NA    интенсив… separ…
#>  3      3 grad1  2       rate o… %/day "\"\""  "\"\""     NA    интенсив… <NA>  
#>  4      4 kdQgr1 150     maximu… %      <NA>   "can be d… NA    максимал… <NA>  
#>  5      5 polmo… 2       the ea… <NA>   <NA>   "can be d… NA    cамый ра… <NA>  
#>  6      6 polmo… 5       the la… <NA>   <NA>   "can be d… NA    самый по… <NA>  
#>  7      7 polko… 15      amount… <NA>   <NA>    <NA>      NA    количест… <NA>  
#>  8      8 polko… 25      amount… days   <NA>    <NA>      NA    количест… <NA>  
#>  9      9 polko… 30      amount… days   <NA>   "can be d… NA    количест… <NA>  
#> 10     10 polgr… 10      mean r… %      <NA>    <NA>      NA    значения… <NA>  
#> # … with 21 more rows, and 1 more variable: role_2 <chr>

Since the number of parameters is large (31), they are organized as list, which is then fed into gr_separate() function. First, you get the params using the gr_get_params() function. The only parameter of the function is `reg``, which indicates the region for which the parameters must be extracted. After you got the parameters, they can be changed by accessing the list elements:

# Расчленение
p = gr_get_params(reg = 'Midplain')
p$nPav = 5
p$prodspada = 85

Next, you can separate the hydrograph:

sep = gr_separate(tab, p)
head(sep)
#>         Date Qin Qbase Quick Qseas Qrain Qthaw Qpb Qtype  Temp Prec
#> 1 1956-01-01 5.2   5.2     0     0     0     0   0     0  -6.5  0.5
#> 2 1956-01-02 5.2   5.2     0     0     0     0   0     0 -11.4  0.8
#> 3 1956-01-03 5.4   5.4     0     0     0     0   0     0 -10.7  0.3
#> 4 1956-01-04 5.4   5.4     0     0     0     0   0     0  -8.1  0.4
#> 5 1956-01-05 5.4   5.4     0     0     0     0   0     0 -11.7  0.1
#> 6 1956-01-06 5.6   5.6     0     0     0     0   0     0 -20.1  0.0

After the hydrograph is separated, it can be summarized in a set of variables:

vars = gr_summarize(sep)
head(vars)
#> # A tibble: 6 x 57
#>    Year Year1 Year2 datestart  datepolend PolProd    Qy  Qmax datemax     Qygr
#>   <int> <int> <dbl> <date>     <date>       <int> <dbl> <dbl> <date>     <dbl>
#> 1  1956  1956  1957 1956-04-08 1956-05-02      24  18.3   467 1956-04-22  7.83
#> 2  1957  1957  1958 1957-03-31 1957-05-01      31  20.3   460 1957-04-08  8.02
#> 3  1958  1958  1959 1958-04-06 1958-05-06      30  27.4   537 1958-04-21  8.28
#> 4  1959  1959  1960 1959-03-31 1959-04-27      27  27.1   406 1959-04-16  7.60
#> 5  1960  1960  1961 1960-03-30 1960-04-26      27  29.5   406 1960-04-15  9.45
#> 6  1961  1961  1962 1961-03-12 1961-04-14      33  18.8   296 1961-04-10  9.81
#> # … with 47 more variables: Qmmsummer <dbl>, monmmsummer <date>, Qmmwin <dbl>,
#> #   nommwin <date>, Q30s <dbl>, date30s1 <date>, date30s2 <date>, Q30w <dbl>,
#> #   date30w1 <date>, date30w2 <date>, Q10s <dbl>, date10s1 <date>,
#> #   date10s2 <date>, Q10w <dbl>, date10w1 <date>, date10w2 <date>, Q5s <dbl>,
#> #   date5s1 <date>, date5s2 <date>, Q5w <dbl>, date5w1 <date>, date5w2 <date>,
#> #   Wy <dbl>, Wgr <dbl>, Wpol2 <dbl>, Wpol1 <dbl>, Wpol3 <dbl>, Wpavs2 <dbl>,
#> #   Wpavs1 <dbl>, Wpavthaw2 <dbl>, Wpavthaw1 <dbl>, WgrS <dbl>, WS <dbl>,
#> #   WgrW <dbl>, WW <dbl>, Qmaxpavs <dbl>, Qmaxpavthaw <dbl>,
#> #   datemaxpavs <date>, datemaxpavthaw <date>, SumProd <int>,
#> #   DaysPavsSum <int>, WinProd <int>, DaysThawWin <int>, CvWin <dbl>,
#> #   CvSum <dbl>, CountPavs <int>, CountThaws <int>

Plot and test

These functions from grwat package allow you to:

Plot separation of hydrograph
Plot interannual changes of key water discharge variables
Plot long-term changes of key water discharge variables
Perform statistical tests on all calculated variables

Graphical functions are based on ggplot2 graphics.

Plot separation of hydrograph

You can plot separations for selected years using gr_plot_sep() function:

gr_plot_sep(sep, 1979) # plot single year
#> Warning: Removed 36 rows containing missing values (position_stack).

gr_plot_sep(sep, c(1994, 2001)) # plot two years sequentially
#> Warning: Removed 24 rows containing missing values (position_stack).

And also multiple years on the same layout:

gr_plot_sep(sep, 2016:2019, # plot four years on the same page
            layout = matrix(c(1,2,3,4), ncol=2, byrow = T))
#> Warning in max.default(structure(numeric(0), class = "Date"), na.rm = FALSE): no
#> non-missing arguments to max; returning -Inf

Interannual change variables

To get the detailed description of available variables you can invoke gr_help_vars():

gr_help_vars()
#> # A tibble: 57 x 19
#>       ID Position Width Source Name   Units Unitsen Readtype Type   Test Desc   
#>    <dbl>    <dbl> <dbl>  <dbl> <chr>  <chr> <chr>   <chr>    <chr> <dbl> <chr>  
#>  1     1        1     7      1 year_… <NA>  <NA>    integer  inte…     0 Номер …
#>  2     2        8    10      1 Year1  Год   Year    integer  inte…     0 Год, к…
#>  3     3       18    10      1 Year2  Год   Year    integer  inte…     0 Год, к…
#>  4     4       28    15      1 dates… Дата  Date    Date     Date      1 Дата н…
#>  5     5       43    15      1 datep… Дата  Date    Date     Date      1 Дата о…
#>  6    57        0     0      0 PolPr… Дней  Days    integer  inte…     1 Продол…
#>  7     6       58    10      1 Qy     м^3/с m^3/s   double   doub…     1 Средни…
#>  8     7       68    10      1 Qmax   м^3/с m^3/s   double   doub…     1 Максим…
#>  9     8       78    15      1 datem… Дата  Date    Date     Date      1 Дата м…
#> 10     9       93    10      1 Qygr   м^3/с m^3/s   double   doub…     1 Средни…
#> # … with 47 more rows, and 8 more variables: Descen <chr>, Group <chr>,
#> #   Winter <dbl>, Chart <chr>, Color <chr>, Order <dbl>, Range <chr>,
#> #   Problems <chr>

Parameters can be statistically tested using test_variables(df, ..., year = NULL, locale='EN') function. Names of the parameters are passed comma-separated in place of .... They are quoted, so you do not need to pass them as character strings, just write their names:

gr_test_vars(vars, Qmax)
#> Warning: `select_()` was deprecated in dplyr 0.7.0.
#> Please use `select()` instead.
#> $ptt
#> $ptt$Qmax
#> 
#>  Pettitt's test for single change-point detection
#> 
#> data:  vl[vl_cmp]
#> U* = 481, p-value = 0.01088
#> alternative hypothesis: two.sided
#> sample estimates:
#> probable change point at time K 
#>                              15 
#> 
#> 
#> 
#> $mkt
#> $mkt$Qmax
#> 
#>  Mann-Kendall trend test
#> 
#> data:  vl[vl_cmp]
#> z = -3.3955, n = 64, p-value = 0.0006851
#> alternative hypothesis: true S is not equal to 0
#> sample estimates:
#>             S          varS           tau 
#>  -587.0000000 29785.0000000    -0.2916775 
#> 
#> 
#> 
#> $tst
#> $tst$Qmax
#> 
#>  Sen's slope
#> 
#> data:  values[fltr]
#> z = -3.3955, n = 64, p-value = 0.0006851
#> alternative hypothesis: true z is not equal to 0
#> 95 percent confidence interval:
#>  -5.418605 -1.560976
#> sample estimates:
#> Sen's slope 
#>   -3.642183 
#> 
#> 
#> 
#> $ts_fit
#> $ts_fit$Qmax
#> 
#> Call:
#> mblm::mblm(formula = eval(frml), dataframe = df.theil[fltr, ], 
#>     repeated = FALSE)
#> 
#> Coefficients:
#> (Intercept)        Year1  
#>    7496.073       -3.642  
#> 
#> 
#> 
#> $tt
#> $tt$Qmax
#> 
#>  Welch Two Sample t-test
#> 
#> data:  d1 and d2
#> t = 3.6606, df = 19.167, p-value = 0.001643
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>   75.52882 276.94261
#> sample estimates:
#> mean of x mean of y 
#>  419.7857  243.5500 
#> 
#> 
#> 
#> $ft
#> $ft$Qmax
#> 
#>  F test to compare two variances
#> 
#> data:  d1 and d2
#> F = 1.2603, num df = 13, denom df = 49, p-value = 0.5375
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#>  0.5776842 3.4623990
#> sample estimates:
#> ratio of variances 
#>           1.260324 
#> 
#> 
#> 
#> $year
#> Qmax 
#> 1970 
#> 
#> $maxval
#> $maxval$Qmax
#> [1] 780
#> 
#> 
#> $fixed_year
#> [1] FALSE
#> 
#> $pvalues
#>   N                                            Variable Change.Year    Trend
#> 1 1 Maximum annual discharge during seasonal flood wave        1970 -3.64218
#>         M1     M2 MeanRatio      sd1     sd2 sdRatio Mann.Kendall Pettitt
#> 1 419.7857 243.55       -42 162.9446 145.144   -10.9      0.00069 0.01088
#>   Student  Fisher
#> 1 0.00164 0.53749

This is an example with three variables selected:

tests = gr_test_vars(vars, Qygr, date10w1, Wpol3)
tests$pvalues
#>   N                                                              Variable
#> 1 1 Annual groundwater discharge ("baseflow") during water-resources year
#> 2 2                   First date of 10-day window discharge during winter
#> 3 3                Seasonal flood runoff (with groundwater and rainwater)
#>   Change.Year    Trend      M1       M2 MeanRatio       sd1       sd2 sdRatio
#> 1        1979  0.11219 7.85396 13.30172      69.4   1.84951   2.32226    25.6
#> 2        1989  0.97421  22-Apr   10-Jun      49.0 133.00000 150.00000    12.8
#> 3        1988 -0.10666 10.2846  6.46593     -37.1   5.44084   3.72020   -31.6
#>   Mann.Kendall Pettitt Student  Fisher
#> 1      0.00000 0.00000 0.00000 0.25596
#> 2      0.47594 0.59462 0.17120 0.49666
#> 3      0.00072 0.00628 0.00206 0.03901

If you want to test all parameters, just skip variable names:

tests = gr_test_vars(vars)
tests$year # this is a change year detected for each variable
#>             Wy          Wpol2            Wgr          Wpol1         Wpavs2 
#>           1978           1974           1977           1988           1977 
#>         Wpavs1      Wpavthaw2      Wpavthaw1           WgrW             WW 
#>           1977           2013           1977           1978           1978 
#>           WgrS             WS             Qy      datestart           Qygr 
#>           1977           1977           1978           1970           1979 
#>     datepolend        PolProd           Qmax        datemax          Wpol3 
#>           1987           1987           1970           1970           1988 
#>         Qmmwin        nommwin      Qmmsummer    monmmsummer           Q30w 
#>           1979           1964           1979           1990           1979 
#>       date30w1           Q30s       date30s1           Q10w       date10w1 
#>           1977           1979           2000           1979           1989 
#>           Q10s       date10s1            Q5w        date5w1            Q5s 
#>           1981           2000           1979           1998           1981 
#>        date5s1    Qmaxpavthaw datemaxpavthaw     CountThaws    DaysThawWin 
#>           2000           2010           1993           1969           1984 
#>       Qmaxpavs    datemaxpavs      CountPavs    DaysPavsSum          CvWin 
#>           1960           1965           1980           1970           1983 
#>        WinProd          CvSum        SumProd 
#>           2000           1991           1988

Long-term changes are tested against breaking year, which is calculated for each variable using Pettitt test. However, if you want to use a fixed year, you should pass the desired breaking year into change_year parameter:

tests = gr_test_vars(vars, Qmax, Qygr, change_year = 1987)
tests$ft # Fisher F tests to compare two variances
#> $Qmax
#> 
#>  F test to compare two variances
#> 
#> data:  d1 and d2
#> F = 1.2603, num df = 13, denom df = 49, p-value = 0.5375
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#>  0.5776842 3.4623990
#> sample estimates:
#> ratio of variances 
#>           1.260324 
#> 
#> 
#> $Qygr
#> 
#>  F test to compare two variances
#> 
#> data:  d1 and d2
#> F = 0.6343, num df = 22, denom df = 40, p-value = 0.256
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#>  0.3117119 1.4015893
#> sample estimates:
#> ratio of variances 
#>          0.6342958

Plot interannual changes

Interannual changes are visualized using gr_plot_vars() function. Its syntax is similar to gr_test_vars() and gr_plot_sep():

gr_plot_vars(vars, Qmax) # plot one selected variable

gr_plot_vars(vars, date10w1, Wpol3) # plot two variables sequentially
#> Warning: Removed 20 rows containing non-finite values (stat_smooth).
#> Warning: Removed 20 rows containing missing values (geom_point).
#> Warning: Removed 2 row(s) containing missing values (geom_path).

#> Warning: Removed 1 rows containing non-finite values (stat_smooth).

gr_plot_vars(vars, Qmax, Qygr, date10w1, Wpol3, # plot four variables in matrix layout
                      layout = matrix(c(1,2,3,4), nrow=2, byrow=TRUE)) 
#> Warning: Removed 20 rows containing non-finite values (stat_smooth).
#> Warning: Removed 20 rows containing missing values (geom_point).
#> Warning: Removed 2 row(s) containing missing values (geom_path).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).

You can add the results of statistical tests to the plot by specifying the result of test_variables() function to the tests parameter. In that case the subtitle with test results will be added, Theil-Sen slope and Pettitt test breaking year are drawn as solid (\(p \leq 0.05\)) or dashed (\(p > 0.05\)) lines:

gr_plot_vars(vars, date10w1, Wpol3, DaysThawWin, Qmaxpavs,
             tests = gr_test_vars(vars, date10w1, Wpol3, DaysThawWin, Qmaxpavs)) # add test information
#> Warning: Removed 20 rows containing non-finite values (stat_smooth).
#> Warning: Removed 20 rows containing missing values (geom_point).
#> Warning: Removed 2 row(s) containing missing values (geom_path).

#> Warning: Removed 1 rows containing non-finite values (stat_smooth).

Note that tests parameter of plot_variables() expects the tests for the same variables as those selected for plotting. If you plot variables A, B, C and supply tests for variables X, Y, Z, they will be added without any warnings, and it is your responsibility to keep them in correspondence with each other

Finally, you can plot all variables by not supplying column names to plot_variables() function. In that case tests (if you want to plot them too) should also be calculated for all variables:

gr_plot_vars(vars, tests = gr_test_vars(vars))

Plot long-term period changes

Long-term changes are the differences between summarized statistics of one variable calculated for two selected periods. Because these statistics reflect the differences in distributions of parameters, grwat visualizes them as box plots using gr_plot_periods() function. The syntax is similar to gr_plot_vars() except that you must provide either tests or year parameter. If both are supplied then tests is prioritized (you can also supply a fixed year when testing variables:

gr_plot_periods(vars, Qy, year = 1978)

gr_plot_periods(vars, Qy, tests = gr_test_vars(vars, Qy))

Multiple plots can be combined on one page using layout parameter:

gr_plot_periods(vars, Qy, Qmax, 
                tests = gr_test_vars(vars, Qy, Qmax),
                layout = matrix(c(1,2)))

To plot long-term changes for all variables just skip variable names in function call:

gr_plot_periods(vars, tests = gr_test_vars(vars))

There is also a small helper function that plots a histogram of minimal discharge month for summer and winter periods:

gr_plot_minmonth(vars, year = 1985)

Reports

To render HTML report just pass separation and variables to gr_report() function, and provide the output file name:

report = paste(getwd(), 'Spas-Zagorye.html', sep = '/')
gr_report(sep, vars, output = report)
browseURL(report)

Spas-Zagorye.html

Hydrograph separation

Interannual changes

Long-term changes

Statistical tests

p-values of statistical criteria
N	Variable	Change.Year	Trend	M1	M2	MeanRatio	sd1	sd2	sdRatio	Mann.Kendall	Pettitt	Student	Fisher
1	Annual runoff	1978	0.07966	22.72886	27.16913	19.5	7.60958	7.52459	-1.1	0.13958	0.14553	0.03142	0.92052
2	Seasonal flood runoff (without groundwater)	1974	-0.10182	10.76172	6.14213	-42.9	5.18134	3.61374	-30.3	NA	0.00421	0.00267	0.05912
3	Annual groundwater runoff	1977	0.13915	8.77161	15.83465	80.5	2.41099	5.01904	108.2	1e-05	0	0	0.00084
4	Seasonal flood runoff (with groundwater)	1988	-0.10188	10.05498	6.41846	-36.2	5.15313	3.59534	-30.2	NA	0.0061	0.00206	0.05039
5	Rain-flood runoff (without groundwater)	1969	0.00062	3.01154	2.86082	-5	4.61235	3.19159	-30.8	NA	0.75438	0.91611	0.07746
6	Rain-flood runoff (with groundwater)	1977	0.02143	5.45264	7.32468	34.3	4.61542	5.16170	11.8	NA	0.223	0.15677	0.61042
7	Thaw-flood runoff (without groundwater)	2013	0.00168	1.11322	2.55362	129.4	1.98440	2.67068	34.6	0.67234	1	0.21075	0.22686
8	Thaw-flood runoff (with groundwater)	1977	0.02143	5.45264	7.32468	34.3	4.61542	5.16170	11.8	NA	0.223	0.15677	0.61042
9	Groundwater runoff during winter	1978	0.0414	2.68291	5.21297	94.3	1.51114	4.64647	207.5	0.00021	5e-05	0.00216	0
10	Winter low flow runoff	1978	0.04513	3.66991	6.61922	80.4	2.96717	6.19911	108.9	0.00218	0.00072	0.01254	0.00062
11	Groundwater runoff during summer	1977	0.08197	5.10734	9.3704	83.5	1.83627	2.94587	60.4	3e-05	0	0	0.02545
12	Summer low flow runoff	1977	0.09643	7.58077	12.10505	59.7	4.72252	5.61700	18.9	0.00532	0.00294	0.00149	0.40664
13	Annual discharge during water-resources year	1978	0.0141	19.64125	23.0792	17.5	6.11827	6.52138	6.6	0.72378	0.39157	0.04257	0.77177
14	First date of a seasonal flood wave	1970	365.22694	01-Apr	28-Mar	-4	7.00000	20.00000	185.7	0.77616	0.59462	0.19563	0.00047
15	Annual groundwater discharge (“baseflow”) during water-resources year	1979	0.11099	7.85396	13.28411	69.1	1.84951	2.32539	25.7	0	0	0	0.25314
16	Last date of a seasonal flood wave	1976	365.06742	03-May	20-Apr	-13	19.00000	24.00000	26.3	0.03894	0.07722	0.0289	0.22037
17	Duration of a seasonal flood wave	1987	-0.2	30.29032	22.84848	-24.6	14.86426	17.08823	15	0.00026	0.00042	0.06739	0.44521
18	Maximum annual discharge during seasonal flood wave	1970	-3.64218	419.78571	243.49	-42	162.94462	145.21444	-10.9	0.00069	0.01088	0.00164	0.53899
19	Date of a maximum annual discharge during seasonal flood wave	1967	365.22727	16-Apr	15-Apr	-1	5.00000	42.00000	740	0.51975	0.55801	0.95763	0
20	Seasonal flood runoff (with groundwater and rainwater)	1988	-0.1061	10.2846	6.4776	-37	5.44084	3.75705	-30.9	9e-04	0.00689	0.00218	0.04425
21	Minimum monthly discharge during winter	1979	0.1239	4.1913	9.52927	127.4	1.98973	2.38341	19.8	0	0	0	0.36874
22	Month of a minimum monthly discharge during winter	1964	0	12-Jun	09-May	-34	54.00000	45.00000	-16.7	0.9386	0.71943	0.13272	0.42549
23	Minimum monthly discharge during summer	1979	0.06406	5.8	9.27561	59.9	5.62737	5.90829	5	0.00115	1e-05	0.02419	0.82616
24	Month of a minimum monthly discharge during summer	1990	0	06-Jun	21-Jun	15	39.00000	44.00000	12.8	0.09557	0.14111	0.16502	0.47698
25	Minimum 30-day window discharge during winter	1979	0.1031	6.19841	10.8255	74.6	2.29125	2.18115	-4.8	NA	0	0	0.76775
26	First date of 30-day window discharge during winter	1977	1.5	17-Apr	14-Jun	58	147.00000	163.00000	10.9	0.0616	0.19193	0.1694	0.64123
27	Minimum 30-day window discharge during summer	1979	0.06228	6.25773	9.75442	55.9	1.70763	2.36964	38.8	NA	4e-05	0	0.11114
28	First date of 30-day window discharge during summer	2000	365.33333	13-Jul	30-Jul	17	29.00000	28.00000	-3.4	0.68392	0.24187	0.03963	0.92991
29	Minimum 10-day window discharge during winter	1979	0.10816	5.27522	10.16122	92.6	1.73743	2.66385	53.3	0	0	0	0.03518
30	First date of 10-day window discharge during winter	1989	0.9881	22-Apr	10-Jun	49	133.00000	150.00000	12.8	0.47594	0.59462	0.17057	0.49787
31	Minimum 10-day window discharge during summer	1981	0.06235	5.93875	8.90868	50	1.58371	1.90763	20.5	NA	3e-05	0	0.34881
32	First date of 10-day window discharge during summer	2000	365.46154	21-Jul	08-Aug	18	26.00000	31.00000	19.2	0.27383	0.07921	0.03672	0.35453
33	Minimum 5-day window discharge during winter	1979	0.11155	4.90696	9.91122	102	1.74827	2.47677	41.7	0	0	0	0.08371
34	First date of 5-day window discharge during winter	1998	0.68819	21-May	30-Mar	-52	146.00000	118.00000	-19.2	0.25602	0.54022	0.12579	0.30373
35	Minimum 5-day window discharge during summer	1981	0.06276	5.81333	8.73474	50.3	1.54071	1.86502	21	NA	2e-05	0	0.33605
36	First date of 5-day window discharge during summer	2000	365.5	24-Jul	11-Aug	18	26.00000	30.00000	15.4	0.21508	0.06023	0.03352	0.38446
37	Maximum thaw-flood discharge	1978	0.03541	21.84435	29.14668	33.4	38.07410	38.38868	0.8	0.54297	1	0.47135	0.99844
38	Date of a maximum thaw-flood discharge	1993	266.27778	30-Aug	04-May	-118	146.00000	133.00000	-8.9	0.0121	0.09772	0.00137	0.63766
39	Number of thaw-flood events	1969	-0.09091	16.07692	11.5098	-28.4	7.46616	6.16887	-17.4	0.01896	0.13894	0.05826	0.33884
40	Number of days with thaw-flood events	1984	-0.53638	92.28571	74.33333	-19.5	29.40755	50.76332	72.6	0.01493	0.008	0.08142	0.00446
41	Maximum rain-flood discharge	1960	-0.16821	125.91617	56.90355	-54.8	98.79488	50.38677	-49	0.49053	1	0.25716	0.02792
42	Date of a maximum rain-flood discharge	1965	365.9099	13-May	14-Jul	62	35.00000	80.00000	128.6	0.25362	0.48353	0.00073	0.0195
43	Number of rain-flood events	1991	-0.05	12.44828	10.2381	-17.8	4.98273	5.02896	0.9	0.24364	0.3658	0.13095	0.94633
44	Number of days with rain-flood events	1970	0.05719	95.78571	107.54	12.3	38.26060	39.61046	3.5	0.81668	1	0.32462	0.94524
45	Relative variation of discharge during winter low flow	1983	-0.00039	0.37243	0.37148	-0.3	0.27190	0.31670	16.5	0.84838	1	0.98974	0.42184
46	Duration of winter low flow period	1968	-0.22967	147.75	130.30769	-11.8	60.62272	68.16432	12.4	0.26575	0.51707	0.39206	0.70438
47	Relative variation of discharge during summer-autumn low flow	1991	-0.0018	0.679	0.57022	-16	0.32442	0.24002	-26	0.35695	0.76859	0.12886	0.10584
48	Duration of a summer-autumn low flow period	1988	0.36549	197	212.8125	8	47.11756	45.36052	-3.7	0.05657	0.26671	0.17637	0.83377

Introduction to grwat R package

Timofey Samsonov

2021-07-09